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We develop and calibrate a new method for estimating the gravitational radiation emitted by 
complex motions of matter sources in the vicinity of black holes. We compute numerically the lin¬ 
earized curvature perturbations induced by matter fields evolving in fixed black hole backgrounds, 
whose evolution we obtain using the equations of relativistic hydrodynamics. The current implemen¬ 
tation of the proposal concerns non-rotating holes and axisymmetric hydrodynamical motions. As 
first applications we study i) dust shells falling onto the black hole isotropically from hnite distance, 
ii) initially spherical layers of material falling onto a moving black hole, and hi) anisotropic collapse 
of shells. We focus on the dependence of the total gravitational wave energy emission on the flow 
parameters, in particular shell thickness, velocity and degree of anisotropy. The gradual excitation 
of the black hole quasi-normal mode frequency by sufficiently compact shells is demonstrated and 
discussed. A new prescription for generating physically reasonable initial data is discussed, along 
with a range of technical issues relevant to numerical relativity. 
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I. INTRODUCTION 

The numerical investigation of relativistic gravity is considered an important complement to both the analytic 
studies of strong field phenomena, and the experimental efforts anticipated to probe such astrophysical regions. With 
the exception of the vacuum two-body problem (i.e., the coalescence of two black holes), all realistic sources of 
gravitational radiation involve matter. The joint integration of the equations of motion for matter and geometry was, 
hence, in the minds of theorists from the very beginning of numerical relativity (e.g., see recollections in ||^). Still, the 
effort focused from early-on on vacuum spacetimes, since those clearly captured, and still do, the least understood part 
of the equation. Today the struggle with the perplexities of the “many fingered time” in relativity may be gradually 
turning in favor of the numerical relativists, at least for certain classes of spacetimes (e.g., see recent work in |^-^). 

In the meantime, largely motivated by observational data on (special) relativistic astrophysical phenomena, the 
numerical investigation of the relativistic equations of hydrodynamics has made signihcant progress. The traditional 
approach to the numerical integration of the hydrodynamical equations started with Wilson’s pioneering work [|j, 
where the equations were written as a set of advection equations. Wilson’s scheme used a combination of upwind finite- 
difference and artificial viscosity techniques to damp spurious oscillations. Despite its non-conservative character, this 
approach has been widely (and successfully) used over the years. Recently, though, high resolution shock-capturing 
numerical schemes (HRSC) have emerged as a solid alternative, which accurately resolves ultra-relativistic flows. On 
the theoretical side this progress was largely based on the skillful exploitation of the hyperbolicity of the system. On 
the practical side, the enormous body of work on numerical non-relativistic hydrodynamics provided a multitude of 
working algorithmic techniques, which were successfully transcribed to the relativistic case (for a review of the current 
status of the field see [||Q|). Clearly, handling the matter degrees of freedom numerically is an increasingly mature 
subject, which will in turn help explore large classes of interesting spacetimes. 

In this paper we start a line of investigation where we will be exploring two and three-dimensional relativistic 
hydrodynamics around black holes, coupled to approximate versions of the equations governing the black hole space- 
time geometry. The aim is to produce simple and robust numerical probes of the gravitational emission of complex 
hydrodynamic motions in the vicinity of black holes. Current and near-future advances in space-borne high-energy 
instrumentation, e.g., the Rossi X-Ray Timing Explorer (RXTE) bring the strong held regions near black holes in 
increased focus. At the same time, some of the most ubiquitous sources of gravitational radiation, anticipated to be 
detected by the Laser Interferometric Gravitational Observatory (LIGO), are expected to be binary systems in which 
a black hole is either present initially, or forms during the merger. The description of the violent merger phase is at 
present understood rather poorly. The combination of theoretical computations with observational output will offer, 
potentially, some thrilling insights into the relativistic regime. To this end, we propose to use special versions of black 
hole perturbation theory, in tandem with full hydrodynamical simulations, in order to estimate the gravitational wave 
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signal from a range of systems. The approach attempts (in effect) to eliminate the need to reproduce and main¬ 
tain numerically the underlying black hole spacetime, complementing the on-going efforts for a complete (non-linear) 
evolution. 

In recent research, vacuum black hole perturbations were integrated numerically as 2-1-1 initial value problems 1^,^. 
This work followed the propagation of generic vacuum initial data through the burst, quasi-normal ringing and power- 
law tail phases of the dynamics, confirming anticipated behavior, even for delicate effects like super-radiance |10f| . The 
concept of numerical integration of linearized black hole perturbations is carried with the proposal presented here one 
step further, by the inclusion of source terms generated by extended objects, in particular from fluid material. 

Materializing the general proposal, we implement, test and use a restricted subset, which still covers new ground. 
More precisely, the hydrodynamical code is restricted to axisymmetry. The fluid configurations we study are pressure¬ 
less gas (dust), albeit evolved with the full hydrodynamical code. The black hole perturbations concern a non-rotating 
hole. We finally focus on the quadrupole part of the radiation. All of the above restrictions can be lifted with straight¬ 
forward extensions of the numerical algorithms. 

The organization of the paper is as follows: in Sec. || we lay out the physical assumptions that underlie our model. 
Next, in Sec. IT we compile the mathematical framework that captures those assumptions. The equations of general 
relativistic hydrodynamics, are briefly reviewed here. We describe perturbation equations with extended sources, 
in particular hydrodynamical, perfect fluid, sources. Specializing the discussion, the inhomogeneous Bardeen-Press 
equation for perfect fluid matter fields is analyzed. Two issues pertaining to the selection of initial data are discussed 
here, first, an approach to minimizing initial radiation content, second, practical strategies for avoiding ambiguities. 
Section IV gives a brief description of the numerical methodology employed for the solution of the hydrodynamical 
equations and the inhomogeneous Bardeen-Press equation. The calibration of the present code is described next. 
The discussion of the results (Sec. ^ begins with a detailed presentation of a typical computation. We then present 
axisymmetric dust computations in a variety of parametric surveys. The conclusions of this work and a discussion 
of technical issues that were raised during this investigation are to be found in Sec. VI. The section, and the paper, 
close with an outline of future directions. 


II. THE PHYSICAL MODEL 

The rapid accretion of large lumps of matter onto a black hole emerges as a fairly generic phenomenon (or at least 
picture) in astrophysics. It occurs when the collapse of a massive star gives birth to a central black hole, as chunks of 
stellar material are being swallowed, while others, prevented by angular momentum, are left orbiting and form a long 
lived torus. Wholesale accretion may also be the end act in a black-hole neutron-star merger, where the neutron star 
might be swallowed whole, or might first disrupt and then accrete in pieces. Or, just the same, in a binary neutron 
star merger, where the final configuration is too massive to support itself against collapse and forms a hole surrounded 
by neutron star matter. 

We outline here the idealizations we adopt in order to get a handle on those intriguing astrophysical events. Those 
assumptions are fairly standard in certain fields, but our numerical approach will be clarified with a recap of the 
main points. In our treatment of accreting systems we take an exact black hole (vacuum) background spacetime as 
the stage on which the accretion play unfolds. The dynamics of an idealized fluid is then explored using the zeroth 
order geometry defined by this spacetime. The stress-energy of the moving fluid generates variable distortions of the 
spacetime, which sufficiently far from the black hole (and the matter source) are interpreted as gravitational waves. 
The core assumption is then that the mass /x, of the accreting fluid, is much smaller than the black hole mass M, i.e., 
AX << M. 

To first order in /x, the matter flow generates metric distortions which feedback as geometric gravitational forces 
acting on the matter. Two of the most important such feedback mechanisms are the self-gravity of the fluid, and 
radiation reaction effects. We will neglect in the current programme both mechanisms, by ignoring the first order 
metric corrections to the fluid equations of motion. The first approximation (i.e., no self-gravity) would be valid 
for motions at a radius R close to the black hole, where tidal forces dominate the fluid self-gravity. Assuming a 
characteristic fluid size r, the radius at which this happens is given approximately hy R k r(M//i)^/^. The second 
approximation (i.e., no radiation reaction) would be valid for sufficiently short intervals of time and weak emission 
of waves. None of those approximations is strictly true in most scenarios of astrophysical interest. Hence we adopt 
them with the understanding that any new results pertaining to realistic systems are subject to scrutiny with respect 
to their actual applicability. Still, it has been a remarkable recent success story, to illustrate the robustness of black 
hole perturbation theory as a quantitative analysis tool 00- 

The remaining dominant effects then, governing the fluid motion in our approximation, are the zeroth-order ac¬ 
celerations and inertial forces encoded in our choice of a static background spacetime, and, additionally the pressure 
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gradients inside the extended object. Neglecting the pressure gradients is also a common assumption in the litera¬ 
ture, which then leads to the “cloud of dust particles” picture and motion along geodesics. Indeed, neither internal 
stresses nor inertial forces will compete successfully against gravity in an advanced stage of accretion, deep inside the 
potential well. In the earlier stages, though, pressure gradients and their interplay with zeroth-order effects may play 
an important role. This is certainly the case for configurations that are in rotating equilibrium, as for example in the 
structure of fluid tori orbiting around black holes. The large number of idealizations we adopted till this point suggests 
that a simple description of the fluid thermodynamics would be sufficient for our purposes. Indeed, the generally fast 
timescale of the collapse reduces our concerns about radiative or dissipative processes. Hence we will work within the 
framework of perfect fluids, and for the most part we will be using in our investigations standard polytropic equations 
of state, even though in the applications considered in this paper we actually assume zero pressure. 

We are still left with important sources of complexity. Whereas the integration of linear equations with potentials 
can be handled in certain cases entirely within analytic techniques fl^ , the source terms from realistic extended 
objects involve steps that appear to be intractable analytically: i) The time evolution of matter fields is generally 
non-linear, with the exception of certain classes of fundamental fields (e.g., scalar and electromagnetic fields). Simple 
initial data for fluids may develop rapidly into solutions with highly non-linear features (e.g., shock formation in 
supernova collapse), ii) The stress-energy tensor is non-linear in the matter fields, iii) The source term generating 
curvature waves depends on rather intricate additions and differentiations of the stress-energy tensor. This reflects 
the dependence of gravitational radiation on higher order moments of the stress-energy. After the buildup is 
complete, those linear operations “separate the wheat from the chaff”. 

Radial flows, i.e., flows with purely radial velocity fields, are a natural subset of the bewildering variety of general 
flows. Their relative simplicity will provide us here with a testing ground for our algorithms and a starting point 
towards more general cases. In the case of dust matter, we may further decompose radial motions into isotropic radial 
flows, where the pattern of fluid velocity does not have any angular dependence (hence all emission comes from the 
accelerated motion of some angularly structured density), and anisotropic radial flows in which the velocity field is 
direction dependent. Finite-sized collections of dust, shaped in the form of stars or shells, falling onto the black hole 
isotropically, have been studied quite extensively ll4| - [l^ . These studies brought to the foreground the fact that for a 
fixed amount of infalling mass the gravitational radiation efficiency is reduced compared to the point particle limit. In 
fact, those studies indicate that the energy emitted by infalling dust will never exceed that of a particle with the same 
mass [ p^ , p^ . The reduction is, intuitively, due to cancellations of the emission from distinct parts of the extended 
object. It appears that the study of increasingly realistic extended objects is crucial for an accurate assessment of the 
waveform of the generated waves in the late stages of a binary merger. 

Within the assumptions mentioned, and the methodological restriction to radial dust flows, our selection of initial 
conditions attempts to cover the range of qualitatively different possibilities. We focus on highly resolved radial 
structures and attempt to keep the angular dependence to the minimal necessary. We anticipate that this trend 
should be reversed in configurations with angular momentum. 


III. ANALYTIC FRAMEWORK 

A. Perfect fluid hydrodynamics in background black hole spacetimes 

The motion of matter fields in a general curved spacetime is governed by the local conservation laws of baryon 
number and energy-momentum 


= 0 , = 0 , ( 1 ) 

where = pu^ is the density current, = phu^u’^ + pg^'^ is the stress-energy tensor for a perfect fluid and 
is the covariant derivative associated with the background metric In those expressions all quantities have their 
usual meaning, that is, p is the rest-mass fluid density, h is the specific enthalpy, defined as/i=l-|-e-|-p//9, eis the 
specific internal energy, p is the pressure and is the fluid 4-velocity. The system of equations must be closed with 
an appropriate equation of state, p = p{p, e). Greek (Latin) indices run from 0 to 3 (1 to 3) and we use natural units 
(G = c = 1) throughout. 

By choosing an appropriate set of variables, the equations of general relativistic hydrodynamics are written as a 
hyperbolic system of balance laws. In ^ this is done by defining quantities which are directly measured by Eulerian 
observers, i.e., the rest-mass density D = pW , the momentum density in the j-direction Sj = phW^Vj and the total 
energy density E = phW"^ —p. In these expressions W stands for the Lorentz factor, which satisfies IF = (1 — 
with , where u® is the 3-velocity of the fluid, defined as u® = vd jW -I- (3'^ ja, where a and /3® are the 
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spacetime lapse function and shift vector, respectively, and 7 ij are the spatial components of the spacetime metric 
where the fluid evolves: 


ds^ = —(a^ — j3ifi'')dt^ + 2j3idx'^dt + 'jijdx^dx^ . 


For a generic spacetime the system of equations we solve reads 

V dt dx'^ ) 


S(w), 


( 2 ) 

( 3 ) 


where g = dei[g^v) is such that 


= aV7, 7 = det(7y), 


and det stands for the determinant of the corresponding matrix. In this work we use the standard form of the black 
hole metric in Boyer-Lindquist coordinates (t,r, d, </>), so the above expressions are specialized accordingly. 

In Eq. (^) the state vector of primitive variables is defined as w = {p,Vi,e) and the vector of evolved variables is 
U(w) = {D,Sj,T), with T = E — D, that is the total energy density subtracting the rest-mass density. In addition, 
the fluxes are 


F‘(w) = (l> (.* - f) 9 

and, finally, the corresponding source terms S(w) are given by 

S(w) = (o.T-- - Cw,) ,o , (5) 

where Ff,^ are the 4-dimensional Christoffel symbols. 

The physical boundary conditions we impose on the hydrodynamic quantities depend largely on the scenario under 
study. In all cases considered in this work we adopt ingoing radial boundary conditions (towards the hole), both at 
the horizon and at the outermost zone. The former condition (at the horizon) is strictly valid as the fluid accretes onto 
the black hole supersonically. The latter condition (outer boundary) implies a continuous inflow of very low density 
matter through the outer boundary of the domain. In the angular direction, axisymmetry considerations prescribe 
the appropriate conditions at the poles (0 = 0 and tt). 



B. Curvature perturbations by hydrodynamical sources 

A direct integration of a perturbative system of the Einstein equations has to follow closely, methodologically, the 
corresponding non-linear approaches, and actually carries a similar computational load, along with gauge and bound¬ 
ary problems. We choose instead to base our approximation of the black hole geometry on the highly developed theory 
of curvature perturbations (the subject is comprehensively exposed in the monumental book of Chandrasekhar p^ ). 
This formalism is teeming with positive features: the evolved quantities are invariant under infinitesimal coordinate 
and tetrad transformations [^ , all possible types of perturbations are neatly captured in the real and imaginary 
parts of the Weyl tensor scalar, and rotating black holes are covered by simple extensions of the non-rotating case. 
Last, but not least, the formalism leads quite naturally to partial differential equations defined in the time domain. 

The dynamics of the Weyl spinor abcd is given generally by 

= 47rV^rc_D)A'B' , (6) 

which shows that spacetime derivatives of the stress-energy tensor act as sources to the tidal field 4/. Those 
evolution equations for the 4* components are valid for any spacetime satisfying the Einstein equations, and reflecting 
their origin in the Bianchi identities, are “already linearized” in the dependent quantities. Their structure is similar 
to the propagation equations for fields of other spin, most notably the Maxwell field. The important difference is that 
the metric, implicit in the connection V;^/, is here part of the theory and must be propagated forward, along with a 
host of other variables in order to close the system of equations. 

The concept of evolving the curvature components becomes radically simpler in the framework of the Newman- 
Penrose formalism, when applied to black hole perturbations ||^,^. For a general class of exact solutions, the 
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equations for the perturbed Weyl tensor components decouple, and when written in the frequency domain even 
separate Here we will use the limiting case for non-rotating black holes, i.e., the inhomogeneous Bardeen-Press 

(BP) equation |^. We adopt the standard form of the equation, i.e., Boyer-Lindquist coordinates and the Kinnersley 
null tetrad (/^, m^). We consider solely the spin —2 field, which encodes, at large r, the outgoing radiation. 

For brevity we use the operators 


d d 

dt~^ i9r* ’ 

(7) 

d d 

dt dr^ 

( 8 ) 

^ d /, * '9 

9 = + s cot 6 » . „ „ , , 

oO sin 0 0 (p 

(9) 

R d /, * '9 

^ = an ^ ^ • fl ’ 

oO sin 0 0 (j) 

( 10 ) 


where 9,9 are the standard eth and eth-bar operators (see e.g., and s is the spin-weight of the field acted 

upon by the angular operators. We use the tortoise coordinate r*, defined by dr^, = {r^/A)dr, where A is the horizon 
function A = — 2Mr with M being the mass of the black hole. We use as dependent variable the complex field 

Y, related to the Weyl tensor tetrad component '£'4 = —C^^pan^m'^nPfrY by F = rt[' 4 . With those choices, the BP 
equation reads 


A 

□F = Stt— r4, 
r 


where the “wave” operator expands as 




The source term of equation (pd]), transcribed from in our notation, is 


1 -- 1 A^ 

T4 = — 7—ySdTnn — - ^Tfnfh 

2r^ 4 

,3r-5M^ ^ , 3r-8M^^ 1 ^ 

“1“ ^9 L — TjYijn -\- /— Q /— L — Uln 

2 r 2 2y'2r^ V2r 


( 11 ) 


( 12 ) 


(13) 


where T.g are the contractions of the stress-energy tensor with the null-tetrad. 

The analytical problem should be clear at this point: given initial data F, F on an initial spacelike hypersurface 
to, and a spacetime function Tn{x^), compute the value of F to the future of to. We further assume that r 4 has 
compact support and we are interested, in particular, in the “time signal” F(t, Tq), where Tq is some remote observer 
location. The limit case where the source is non-zero only along a timelike trajectory has been studied extensively with 
semi-analytic methods, mostly in connection with estimates of inspiraling compact binary waveforms (see e.g., ^,^^). 

The higher derivatives of metric quantities involved in the curvature approach to perturbations limit the amount of 
geometric information that can be readily extracted in time-domain integrations of the decoupled equations. Assuming 
the components of the Weyl curvature tensor have been reconstructed, the complete metric information is in principle 
extractable |p0| but in practice may be cumbersome to compute, especially in the rotating black hole case. This is, 
of course, not a problem at infinity, where the gravitational wave signal is to be extracted. 

The asymptotic behavior of the field F along ingoing and outgoing null directions is important for numerical 
purposes, as it reflects the radial behavior of generic solutions starting with data of compact support. The asymptotics 
of outward propagating solutions to the homogeneous equation are given by 


lim F ~ 1, lim F ~ 1, 

r*—»-+cxD r *—»■ —00 


while inward propagating solutions behave as 


lim Y ^ r'^, lim F 

r *—>-+oo r*—* — oo 


A2 . 


(14) 


(15) 


The exponential radial decay (in r*) of inward propagating solutions near the horizon implies F = 0 there. In 
practice, assuming pure ingoing waves, (T_F = 0), at a finite but very small distance from the horizon works well. 
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The asymptotic dependence is not as favorable in the outer boundary (this issue is discussed further in Sec. VI). In 
practice we consider only the domain of dependence of the initial data. In two-dimensional simulations the symmetry 
axis is handled with simple regularity conditions. 

The analytic formulation of the evolution equation for the s = 2 component of the Weyl tensor is entirely analogous 
to the s = —2 case and it would appear that similar algorithms would be effective. Nevertheless, the limiting behavior 
of this fields near the horizon leads to numerical problems (already present in radial integrations in the frequency 
domain |2^). 

The contractions of the perfect fluid stress-energy tensor with the null tetrad read r„„ = phu^, Tnm = phunU and 
Tfhfh = phu^, where we introduced the spin-weighted functions = u^n^,u = and u = u^fh^ with spin- 

weight 0,1,-1 respectively. Some interesting simplifications that have emerged from the use of the Newman-Penrose 
formalism are notable here. The entire physical content of the flow, to the extend that it affects 'I' 4 , is captured by 
the four scalar functions (p, (u is complex), where p = ph. 

We single out the sub-case ^ 0 ,m = 0, which corresponds to purely radial motion. For such motions only the 
Tnn term in Eq. dO) remains. It is worth noting that the other extreme geometry, of purely rotational motion^ as 
is for example approximately the case in a neutron-star black-hole binary, does not zero-out the Tnn term, as Tnn 
is bounded below by the rest mass density. In such rapidly rotating material the emission will be governed by the 
combined effect of the time variability in the angular structure of T„„ and the strong fluid currents encoded in Tmm- 

In the case of purely radial motion, we have the further subdivision according to whether the radial velocity 
{v = \/^ijV’-v^) is isotropic or not, i.e., whether = 0. In the general (anisotropic) case the radial motion source 
term expands as 


SSTrin = 99/5 + 4u„9p9M„ -I- 2 m„/599u„ -I- 2p9M„9u„. (16) 

The contributions to the source term can now be analyzed with respect to their angular geometry. In the isotropic case 
all wave emission is induced by the angular profile of the density (99/5), although its intensity will also depend on the 
magnitude of the radial velocity (the factor). It is easily seen that in this case, the source separates completely, and 
different multipoles of the density directly excite the different multipoles of radiation. This can be shown explicitly 
by e.g., raising the spin of the equation twice and then expanding the function R = 99E in spherical harmonics to 
reach the form 


^iRlm — StT — Un [l{l ~ 1)(^ + 1)(^ + 2)] Pin 
r 


(17) 


where 


^ ^ 4(r-3M)^ 6MA A,,, 

□/ — L_L+---L+ -I-r- jl{l + 1). 


(18) 


Assuming angular dependence for the radial velocity generates angular dependence for and hence activates all 
terms in Eq. (|l6|). For example, a spherically symmetric shell would still generate a non-zero source term if the radial 
velocity of the fluid had at least a dipole component. Expanding the expression of the fourth term in Eq. ( 0 ) in 
spin-weighted harmonics shows that this combination produces a spin (-2) quadrupole. Similarly, a dipole density 
distribution is seen to couple with a dipole moment in the velocity to produce an additional lowest order contribution 
to the quadrupole emission. This interesting non-linear complication, introduced by anisotropy, does not appear to 
have been analyzed so far in studies of matter collapse onto black holes. The effect of such leading order anisotropies 
(i.e., dipole) will be probed in section One final comment on the structure of the source term for fluids, in particular 
radial flows, concerns the types of perturbations excited by the sources. We see that radial motions in axisymmetry 
generate only real perturbations, the imaginary part of the field being associated with non-axisymmetric patterns of 
density/radial velocity. 


C. Initial data for the perturbation fields 

We outline here our procedure for selecting initial data for the perturbation fields. The issue can be of outmost 
relevance for numerical investigations, as spurious behavior may completely mask the sought-after signal. It will be 
demonstrated in Sec. ^ that our concrete computations are performed in a way that bypasses the issue. Here we 
summarize and discuss the evidence we accumulated in the process of dealing with this issue. 

Although 'I '4 is typically associated with propagating gravitational radiation, it is generally non-zero and non¬ 
propagating in regions occupied by matter. For example, in equilibrium matter configurations around black holes. 
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which are possible with support from angular momentum, say in the form of tori, the BP equation admits also 
static solutions. Physically, those solutions represent the first order deviation of the curvature of the equilibrium 
system from that of the underlying vacuum spacetime. Clearly, in this case, the initial data selection procedure 
should accommodate for this static, “Coulombic”, part of ' 1 > 4 . In this paper we focus exclusively on purely radial 
motions with zero angular momentum. Here, there can be no true static solutions. Still, the initial distribution 
of density should be complemented physically with a non-zero curvature, which, without being static, should only 
change appreciably in timescales associated with those of the flow (which we can arrange initially to be very slow). 
Prescribing values for the initial curvature Y, Y, once an initial distribution of T 4 is given, such that there is no 
“unseemly” initial burst of radiation, constitutes the problem of specifying “physical” initial data. 

We investigate the possibility of utilizing the concept of an “approximately” static curvature for setting initial data. 
We computed such solutions by considering the elliptic problem obtained from the BP equation, upon assuming static 
curvature fields. In regular integrations we achieved an entirely equivalent effect by considering the field Z = Y and 
the corresponding equation 


□ Z = Stt—T 4 , (19) 

r 

which is straightforward to obtain from the time derivative of Eq. & above, since all operators commute with dt- 
In all the runs shown below, the initial wave data are assumed to be zero throughout the domain, i.e., Z = 0, Z = 0. 
When using such data, the development of non-zero values for the time derivative Z of the tidal field Y is now 
associated with the time derivative of the source term T 4 . Since those time derivatives can be made zero initially, 
one would expect intuitively that this prescription would better control the response of the curvature to the slowly 
accelerating motion of the matter. We find indeed that this approach significantly reduces the unphysical initial 
radiation content of the data compared to setting zero data on Y. 

The radial infall scenarios are generally weak emitters, and hence in the configurations investigated in this work the 
likelihood of contamination is high. Indeed we find that even with the significant reduction achieved with the above 
recipe, several configurations produce real signals with amplitude below that of the initial burst. In the computations 
below, we opt to leave in the past any ambiguities associated with the choice of initial data. We run configurations 
which, by construction, do not emit appreciably except at late times, hence isolating the two problems. A recent 
study 1 ^ illustrates another clear case of significant radiation hidden in the choice of initial data. 


IV. NUMERICAL METHOD 


A. Integration of the hydrodynamic equations 


We solve the equations of relativistic hydro^namics on a discrete numerical grid with a state-of-the-art HRSC 
scheme. The procedure has been described in Q and successfully applied in We briefly summarize pertinent 

details in this section, referring the interested reader to for further references, details on code calibration, numerical 
tests and algorithmic issues. 

Assuming axisymmetry (2D) and Boyer-Lindquist coordinates (t, r, 6, (p), the vector of evolved quantities, is updated 
from time level to according to the following conservative algorithm 


U»+i ^ u” ■ - — 

Ar 




At 

AO 


■ hj+i 


1 -F 




AfS, 


( 20 ) 


where At = — t", and Ar and A9 indicate the radial and angular grid spacing, respectively. In addition, indices i 

and j label the radial and angular zones, respectively. In Eq. j^), and are the mean values of the state and 
source vector in the corresponding two-dimensional cell, while ^ and F^ are the numerical fluxes which are 

computed at the cell interfaces. These fluxes are calculated using an approximate (linearized) Riemann solver built 
upon the characteristic information of the system Q . The canonical flux-formula we use reads 

Fi±l/2 = 2 + ^i±l/2(wL) — ^ |A(„) |Au;(„)r(„)^j±l/2^ , (21) 


where wl and wr are the values of the primitive variables at the left and ri ght sides of a given cell interface. They 
are derived by a monotonic linear reconstruction of their cell-centered values [p3[. This procedure ensures, in absence 
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of shocks, second-order accuracy in space. In Eq. (pl|), {A„, are the eigenvalues and right-eigenvectors of 

the Jacobian matrices of system (||), evaluated at the cell-interface using the arithmetic mean of the left and right 
states. The jumps of the characteristic variables across each characteristic field, {ALUn}n=i,.., 4 , are given by 

4 

U(wr) - U(wl) = AS3„r„. (22) 

n —1 


The code makes use of a third-order TVD (total-variation-diminishing) Runge-Kutta scheme to perform the time 
integration algorithm of Eq. (p0|), which is done simultaneously in both spatial directions. Finally, a one-dimensional 
Newton-Raphson iteration leads, at each time step, from the evolved variables to the primitive ones 

The well known coordinate singularities of the Boyer-Lindquist coordinates introduce some awkwardness in numer¬ 
ical work. Pure inflow conditions for the fluid cannot be applied arbitrarily close to the horizon, as the coordinate 
velocities become singular there. In the scenarios considered here, the supersonic nature of the flow ensures that 
simple inflow conditions would work well. We have tested that posing the boundary conditions closer to the horizon 
(a range from r* = —2 to = —3) does not affect the flow pattern and the corresponding waveforms. This issue is 
discussed further in Sec. VI. 


B. A simple differencing scheme for the Bardeen-Press equation 


Successful finite difference discretizations of the equation in the homogeneous case have been published in [|j. 
We present here a somewhat simpler algorithm which we implemented and tested in two different guises: i) A 2-1-1 
finite difference computation in which angular derivatives are computed numerically, ii) A multipole decomposition 
approach leading to systems of (decoupled in the non-rotating case) U-l equations. Whereas the first approach easily 
generalizes to rotating holes and higher dimensions, efficiency concerns led us to focus here on the second approach, 
i.e., separating into multipoles, and considering only the low-lying ones. The natural extension of the multipole 
decomposition to the rotating black hole case is through a pseudo-spectral algorithm for the angular terms. 

The BP equation, Eq. (0). is discretized on a uniform tortoise radial grid using a straightforward implementation 
of the three-level leapfrog method. The radial and angular derivatives are approximated as 


[Yoo]l^ = - 2Yi:^ + , 

[YAl, = - E.",_i))/(2A0), 

[Yr,]lj = i,,)/(2Ar*), 

lYr,r,]Y = (^.+10 - i.,)/(Ar.)". 

The time derivatives are approximated as 

[Yu]Y = (1)71 - 2yA + y,7i)/(At)2, 
[Y,t]Y=iK;^-Yl-^))/{2At). 


(23) 


(24) 


Substituting those approximated expressions in the BP equation and solving for the discrete variable k)7^ gives an 
explicit update routine, which is stable subject to the usual CEL condition. 

A direct implementation of the scheme ( p^ , p^ leads at late times (of the order of hundreds of M) to a high 
frequency instability which is seen to emerge from the area around r = 3M. Upon performing a frozen coefficient 
stability analysis of the BP equation, in this region, it can be seen that local propagating modes are unstable even 
in the continuum limit. This instability is inherited to the discrete equations, where, in contrast to the analytic 
canceling of the blowup as local modes propagate away from the region, the numerical modes grow uncontrollably. 
The localized and high frequency nature of the instability leads to simple cures, e.g., by using a numerical scheme 
with higher intrinsic viscosity, or by explicitly adding some artificial viscosity to the discretized equation (the option 
we followed here). 

The inner boundary conditions are those of pure ingoing waves. If those conditions are applied sufficiently close to 
the horizon, this approximation of Eq. is very accurate. In terms of the Boyer-Lindquist r coordinate our typical 
location of the inner boundary corresponds to a coordinate distance of |r — r_|_| « 10“^^, leading to A^ « 10“^®. 

For reasons of numerical efficiency, we have also developed a finite difference implementation of the BP equation 
which explicitly separates Y into spin harmonics. In the non-rotating black hole case each multipole evolves inde¬ 
pendently, under the influence of the corresponding multipole source term, which we construct numerically from the 
function T 4 by performing a projection onto the corresponding spin harmonic: 



( 25 ) 


Tim{t, r*) = j dfl _2YimT4:{t, r^,9, </>). 

Reflecting our physical assumptions and the mathematical model, the coupling of the numerical algorithms is one- 
directional, i.e., numerical information generated by the hydrodynamical code is converted into appropriate numerical 
sources to the BP equation. The generated waveform does not feedback into the fluid motion. An important feature 
of the numerical coupling is that the numerical evolution of the two systems is not done on grids of the same radial 
extend. As already mentioned, the boundary conditions for the wave equation are those of pure ingoing waves, and 
can be applied (and work well) sufficiently inside the potential well. A value of r* = —10 is more than sufficient and 
we typically use r* = —50. In contrast, the inner boundary for the fluid must be kept at somewhat larger radius. 
The situation is reversed in the outer boundary: Appropriate treatment of the asymptotic region is much easier for 
the hydrodynamic quantities than for the conformal held Y. This leads to the adoption of a grid considerably more 
extended in radius for the wave compared to the hydro grid which is kept to the minimal sufficient range. Finally, 
the computation of gravitational wave energy E is performed with the integration of the outgoing signal in time at a 
fixed observer radius, using a second order accurate summation over time and angles, and we use the baryonic mass, 
integrated over the initial hypersurface, as an estimate for the total mass /r. 


C. Calibration of the coupled algorithm 

We demonstrate here convergence results of the gravitational waveform for increasingly finer grid resolution. Al¬ 
though a formal convergence analysis can be done using any finite integration time, even if small, we choose instead 
to check convergence for our typical total evolution times, i.e., of the order of half thousand M. This raises the 
computational requirements for performing the test, but we considered it necessary for properly assessing the validity 
of our results, for the following reason: the timescale of variation of the solution is initially very slow and gradually 
accelerates as the material falls into the hole. A convergence analysis for a small total time will not cover quantita¬ 
tively the most demanding part of the computation, which may span a total time of only a couple of decades in M, 
occurring several hundreds of M after the computation commences. 

Figure ^calibrates the resolution requirements of the wave code in propagating radiation on a black hole background. 
What is shown is the percent relative error e in the total outgoing flux as a function of grid size, i.e., e = {En — Eoo)/Eao, 
where En is the emitted energy for a given radial grid size n and Eoo is the energy out-flux estimated using the finest 
grid. The flux is integrated at r* = lOOM up to a final time t = 400M. The “vacuum” initial data, for this set of 
runs were chosen as 

Z = e-°-2(”*-25)%in2 0, (26) 

Z = 0. (27) 

The resolution requirements are modest, the slope of convergence of the total energy is at about 2.2. The convergence 
curve overestimates the rate when the grid size becomes comparable to the reference grid, so for the computation of 
a numerical convergence rate we dropped the second finest resolution. 

In Fig. 1^, the considerably more demanding convergence behavior for a complete accretion/emission event is demon¬ 
strated. We consider a typical Gaussian quadrupolar shell, fixing the relevant parameters to an initially stationary 
configuration located at r* = 35M, and width k = 2. Details about the initial data and the waveform generated by 
the infalling matter are given in Sec. ^ As with the free wave propagation, we measure the energy flux at r* = lOOM 
up to a final time t = 400M and compute the total energy. The wave initial data are zero for both the held and its 
derivative. The radial extend of the hydrodynamical grid is between —3 < r* < 50. The wave grid extends between 
-50 < r* < 500. 

Resolutions of 0.2M, while more than adequate for capturing typical propagating waves and their interaction with 
the black hole potential, seem to be rather inaccurate in resolving the peak emission event. The plot suggests that 
in terms of adaptive grids (with a refinement ratio of two), (e.g., P^) about two levels of refinement around the hole 
would substantially decrease the computational demands of three-dimensional calculations (assuming that a uni-grid 
resolution of 0.2M can be achieved far from the hole). 


V. RESULTS OF AXISYMMETRIC RADIAL INFALL SIMULATIONS 
A. Radially Imploding Shells: Isotropic velocities 
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1. A typical evolution 


The different initial data sets we explore are parametrized by: i) the location of the shell, vq, ii) its peak density, 
Pmax, iii) its peak velocity Vmax and iv) its width parameter, k, 

P = + (28) 

V = -Vmaxf{0) , (29) 

where po is a given background density, of very small value (typically three to four orders of magnitude smaller than 
the peak density). The angular structure encoded into the functions g{6), f{9) is modelled using simple spherical 
harmonics of low multipole number, e.g., ^ = 0,1, 2. As mentioned before, we restrict our discussion to dust {p = 0). 
The sign of the velocity reflects our further restriction to inflows. 

We present first an overview of a typical evolution. The initial distribution of matter density is described by a 
shell with quadrupolar angular shape, i.e., g{0) = sin^ 9 (see the upper left panel in Fig. |^). The initial velocity Vmax 
is assumed zero. Under the action of the black hole potential, the flow will develop into an isotropic radial velocity 
field. The quadrupolar angular dependence for p is the simplest one that would lead to emission of radiation under 
isotropic radial infall conditions. The shell is initially centered at ro = r* = 35M. The Gaussian parameter is A: = 2. 

The upper panel of Fig. ^ shows the wave signal, more precisely the quadrupole part of Z normalized with the total 
mass of the shell, reaching a remote observer (located at r* = lOOM) as a function of coordinate time. The existence 
of two bursts is the first striking observation. The second burst is more than two times stronger, but the wavelength 
appears to be comparable. To bring out the full range of features, we show in the lower panel the logarithm of the 
same signal. The nature of the signal is now clearer: besides initial low frequency features ahead of each burst, both 
bursts are mostly in the quadrupole ringing frequency of the black hole (17M). 

The existence of two bursts asks for some clarification. To start, we notice that 140M is approximately the time it 
takes for initial radiation to propagate from r* = 35M to about r* = 5M, scatter off the combined angular-momentum 
black hole potential, and reach back to the observer at r* = lOOM. The second burst occurs much latter. Starting 
from rest at r* = 35M, the free fall time of a particle would be around 200M, which combined with the lOOM to 
reach the observer, gives us the arrival time of the interesting emission, associated with the shell crossing the peak of 
the gravitational potential and subsequently accreting onto the hole. 

A spacetime diagram (Fig. |^) of the whole evolutionary sequence makes the previous discussion fairly transparent. 
The use of the tortoise coordinate implies that light rays in the diagram are (modulo graphic distortions) propagating 
at 45° angles. Some guidance to the plot may be worthwhile: the radial coordinate is truncated at the fairly small 
value = 50M to bring out the interesting features. The field variables are sampled at the equator {9 = tt/2). We 
show contours of equal amplitude. In the left panel the quantity contoured is the logarithm of Z, whereas the right 
panel depicts the logarithm of the density. The right panel illustrates the infall of the shell, imperceptibly at first, then 
increasingly faster. Notice the wide spread (in coordinates!) of the accretion event. The left panel shows the response 
of the black hole curvature to the trajectory of the shell. There are three contour levels illustrated, at logarithmic 
amplitudes of —6,0, 2.5 respectively. The —6 level is seen to form a light cone emanating from the initial location of 
the shell. We can clearly see the cone of influence reaching inwards to the horizon region, and exciting quasi-normal 
mode (QNM) ringing. The features of the ringing are traced, e.g., around time t = TOM, by the 0 contour level. As 
this initial burst subsides exponentially, the shape of the infalling shell is emerging, through the imprint it leaves on 
the curvature (still 0 level contours, starting at around t = 120M). At t = 210M the shell is rapidly accreting, its 
quasi-static curvature pattern giving way to emission, which is again seen to be predominantly in QNM form. 

The important fact we deduce from this figure is that the emission we observe is independent of the choice of initial 
wave data provided the collision occurs after a total time U + U, where U is the light crossing time of the combined 
object and U is the time interval in which the black hole ringing decays significantly below the level of the subsequent 
emission. With a conservative assumption of similar amplitudes for unphysical and physical emission, the exponential 
decay of the QNM suggests (see e.g., lower panel of Fig. that U = 50M should give two orders of magnitude 
leeway. 

In Fig. ^ we plot the spatial density distribution at four different times. The symmetry axis 0 = 0 is horizontal, and 
the location of the black hole is depicted with the small black disk at the center. For the plot we use a simple mapping 
of a meridian plane of the Schwarzschild black hole into “Cartesian” coordinates: x = r cos0, y = rsm9. The initial, 
radially very compact, profile of the density is shown at t = 0, followed by the gradual acceleration and infall shown 
at times t = lOOM and t = 200M. The final state is a relic (non-radiative) spherical pattern, corresponding to the 
low background density that “fills the vacuum” surrounding the shell. The t = 200M panel is a four-fold zoom into 
the interesting region. 

Similarly, we plot spatial patterns for the response of the curvature in Figures and |^. The initial burst and 
associated ringing are shown at t = 50M and t = lOOM. At t = 150 those features have subsided sufficiently to reveal 
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the quasi-static field surrounding the shell and being dragged with it towards the hole. At t = 200M, the density 
is near the peak of the potential and the emission is reaching its climax. At t = 250M the excitement is already 
gone: the remaining integration through t = 400M, only obtains the slow decay of the black hole quasi-normal mode 
ringing. 


S. Dependence on shell width, velocity and position 

The narrow radial profile shown in Fig. was picked to produce a loud and clear bang. We illustrate here the nature 
of the emission as a function of radial width, by varying the parameter k of the density data. The initial physical size 
of the shell, i.e., its thickness in terms of proper radial distance is given approximately by L* = 2ac/Vk, where ac is 
the lapse at the center of the shell. This length defines the location where the shell density falls to the 1/e of its peak 
value Pmax- All other parameters stay the same. 

We perform a number of runs, with k parameters ranging from 0.1 to 100 and compute f^ ] the total emitted energy 
by integrating in time the outgoing flux at the location of the observer (r* = lOOM). The dramatic rise in efficiency 
with increased compactness is evident in Fig. and spans several orders of magnitude. In the extended shell limit, 
the energy is shutting off, apparently as a power law of L*. The slope of the dotted line in the insert is very close 
to —2.4. In the infinitesimally thin limit, the energy asymptotes to a finite value, about a third of the point particle 
upper limit. This is consistent with semi-analytic results for oblate and prolate spheroids in |Q. From the insert we 
deduce that there is a characteristic, initial, radial size, of about L* = 1, beyond which no significant efficiency gains 
are made. 

The velocity of the infalling matter, at the time of crossing the potential peak, may be expected to be another 
factor governing the efficiency of the emission, intuitively anticipating larger impact velocities to radiate more. In 
Fig. I we highlight this aspect by plotting the maximum signal amplitude as a function of initial velocity. This is 
done for initial velocities up to 0.1 times the speed of light. The amplitude for the largest initial velocity is more than 
twice that for starting at rest. Fitting a law to the velocity dependence seems problematic without a more delicate 
analysis. The stumbling block is the fact that starting from rest at a finite distance already corresponds to large 
velocities at impact, hence the offset and slope at u = 0 of our diagram are fairly arbitrary and do not capture the 
asymptotic behavior of small velocity. The same figure displays a parallel parameter survey for shells falling from a 
larger distance (r* = 40M). The Gaussian parameter k has been kept the same. This explains the, at first confusing, 
weaker emission from the more remote shell location. Clearly, compactness is overpowering velocity (at least for the 
range of configurations probed here) in controlling the efficiency. 


3. Excitations of QNM emission 

Fig. ^ illustrates the qualitatively different response of the black hole to the infalling shell, depending on the 
shell width. The left column depicts the waveform for four select values of k (increasingly thinner shells from top 
to bottom). The large increase of per unit mass emission is evident. The apparent wavelength of the emission is 
also seen to be increasingly smaller. The right column takes the natural logarithm of the signal, to bring out the 
qualitative changes brought about by the changing shell profile. With progressively more compact shells we see first 
a marginal (second panel), and then a more clearcut (third panel), excitation of the QNM ringing [^. In the last 
panel the ringing is dominant from t = 315M onwards. Hence, as intuitively anticipated, a broad initial distribution, 
emits less, in longer wavelengths, and correspondingly does not excite the lowest QNM frequency of the black hole as 
much. 

The importance of QNM waveforms in the interactions of black hole with external matter was illustrated early 
on More detailed numerical work corroborated the case for QNM excitation in stellar collapse. On the other 
hand, calculations based on point particles scattering of a black hole ig,il showed excitations induced by an external 
object that do not produce significant QNM ringing. In subsequent work the amount of excitation of the QNM 
modes during core collapse has been examined using both analytic model problems and further simplified numerical 
calculations in the spirit of Q . 

Despite the considerably different setup, it appears that our findings support the qualitative and intuitively agreeable 
statement, made in that QNM excitation is induced by the curvature profiles that have spatial wavelengths 
comparable to the width of the black hole potential. In our setup, the width of the shell controls the shape of 
the quasi-static curvature profile. Modulating the profile using the width parameter, strongly affects the degree of 
excitation of the QNM signal. 
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B. Radially Infalling Shells: Anisotropic Velocities 


One of the physical objectives in this paper is to approximately map the role of the fluid source term in the case 
of radial infall. It seems that anisotropic velocity profiles have not been investigated in the literature to the extent 
that isotropic flows have. From the structure of the source term we see that in this case the source non-linearity of 
the right-hand-side of the BP equation becomes important. Quadrupolar (and of course higher) contributions may 
emerge from couplings of angular structure between velocity and density profiles. 


1. Spherical shell collapsing onto a moving black hole 

The geometrically simplest, dipole, anisotropy in the velocity field occurs naturally for matter motion with net linear 
momentum with respect to the black hole. One interesting scenario is the accretion of the layers of a star onto a core 
that has acquired net momentum (kick velocity) during the collapse. To simulate such a configuration we consider a 
spherical Gaussian distribution of density centered at r* = 35M (simulating a dense layer of material, with its various 
parts initially at rest with respect to each other). We assume that earlier events have given the core of the star (now 
collapsed into a black hole), a linear velocity of about 0.1c. In practice, the computation is done in the rest frame of 
the black hole, and the linear momentum of the fluid translates to a radial velocity profile f{6) = cos0. The value 
we adopt for the linear velocity is generous compared to the typical observed kick velocities of neutron stars, which 
reach, at best, one-thirtieth of this speed. The unphysically large value is adopted here for illustration purposes. 

Fig.0 shows the development of the density at intervals of lOOM. The hole rushes up the symmetry axis, towards 
the positive hemisphere (right side of the panels). By t = lOOM the shell has been distorted appreciably, as parts of 
it experience increased acceleration due to the proximity to the approaching hole. By t = 200M, the hole has gone 
right through the shell, and we can see concentrated accretion to occur at the equatorial plane. At this time most of 
the gravitational wave emission has occurred. In the final snapshot (at t = 300M) we can see the slow accretion of 
the remaining shell material from the back side of the black hole. 

The time profile of the quadrupole signal component is essentially identical with the high compactness shells we 
considered before (see e.g., last panel of Fig. The energy released in this event as Z = 2 radiation is measured at 
2.9 X The numerical error, estimated from the radial calibration plot (Fig. ||), indicates that the numerical 

energy estimate is accurate to about 2%. 


2. Pancake collapse 

In the case of rotational collapse one would expect more involved anisotropic velocity profiles. We do not simulate in 
this work rotational configurations, but will attempt to mimic a differential fall using an anisotropic velocity pattern. 
We hence assume again a spherical Gaussian distribution of density centered at r* = 35M and take f{6) = sin 6*. 

Fig. 1^ shows the development of the density at intervals of 50M. The shell collapses onto the black hole with 
the equatorial region falling in first, while the polar regions only gradually accelerate. By t = 150M the equatorial 
zones have accreted and further material keeps falling in from larger latitudes. The energy released here (again in the 
quadrupole mode) is measured at 2.4 x 10“^/r^/M. The accuracy is again at the 2% level. 

Gomparing the efficiency figures for the two anisotropic scenarios considered here we notice the order of magnitude 
difference. Given that the shell material has for both cases similar peak radial velocities on impact, it is likely that 
the key factor in determining the efficiency is what fraction of the total mass accretes with velocities close to the peak 
and radial profile sufficiently compact. 


VI. DISCUSSION 

We developed a new method that lies in the borderline between fully non-linear numerical integrations of the 
Einstein equations coupled to matter fields, and the test particle semi-analytical investigations. Our primary purpose 
is estimating the gravitational radiation emitted by complex motions of matter sources in the vicinity of black holes. 
We opted to approximate the first order deviations from the exact black hole geometry using curvature perturbations 
induced by perfect fluids whose non-linear evolution is integrated using a hydrodynamical code. We used typical infall 
computations to demonstrate the second order convergence of the signals and calibrate the energy estimates over a 
very large amount of evolution time. 
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We studied dust shells falling isotropically from finite distances, with variable shell shapes and initial velocities. We 
explored the dependence of the energy emission on the shell thickness and initial velocity, fitting the numerical data 
when possible. We demonstrated the gradual excitation of the black hole quasi-normal mode frequency by sufficiently 
compact shells, confirming ideas about the mechanisms responsible for producing strong QNM signals. We extend the 
existing literature on radial infall of dust shells by exploring the case of anisotropic radial infall and the dependence 
of the signal on the degree of anisotropy. The anisotropic scenarios we investigate include the collapse of an initially 
spherical layer of material onto a moving black hole and the anisotropic “pancake” collapse of a spherical layer. In 
broad terms, in the case of radial infall onto a non-rotating black hole, the dominant factor controlling the strength 
of the emission is, quite generically, the radial compactness of the source. 

We introduced and explored a prescription for minimizing the initial radiation content which appears to perform 
rather well in the current context. The prescription concerns the selection of the freely specifiable data, on the 
initial hypersurface, and essentially consists of solving the static problem for the electric and magnetic parts of the 
curvature tensor, derived from their corresponding evolution equations. In the presence of non-zero matter sources, 
the procedure generates the static electric and magnetic curvatures associated with density distributions and fluid 
currents. The approach can be transfered, at least in principle, to systems of non-linear evolution equations that 
incorporate some, or all, components of the Weyl tensor in the evolution variables. 

Despite the reduction in the amplitude of the initial burst, the very low efficiency of the sources we simulated 
in the present work called for more drastic action. The remnant radiation in the initial data has been dealt here 
in an unambiguous, brute force, way: the computation simply marches along for sufficient time, past any shadows 
cast by the initial data prescription. This illuminates an important facet of numerical relativity computations: In 
the absence of convincing arguments demonstrating minimal contamination from the initial data, trustworthy signals 
can only be obtained by evolving for sufficiently long times. This allows the first burst of radiation to decay and 
reach negligible levels which do not obscure the real signal. If this strategy is adopted, a reduction of the amount of 
radiation in the initial data has positive practical consequences, as it reduces the computational resources required 
for unambiguous identification of the signal, especially in three-dimensional (3D) simulations. It is reassuring though, 
that the problem appears to become less relevant in precisely such 3D long term integrations of strong emitters, such 
as particles orbiting around black holes |^]. In that case, the strong quadrupolar emission of the binary appears to 
be quickly dominating the initial burst. 

The flow of matter in the vicinity of the black hole, for a major initial part of the computation, involves timescales at 
least one order of magnitude slower than the timescale of propagating gravitational waves. This circumstance may be 
encountered quite frequently when attempting to simulate realistic dynamical events, as the progenitor (slow) system 
must be modelled for a sufficient amount of time to ensure proper initial configuration. This disparity in timescales 
implies that explicit methods of integration of the geometry (as the one employed here) are constraining the efficiency 
of the computation. Implicit integrations of the wave degrees of freedom may offer a more economical alternative, by 
allowing an integration at the hydrodynamical timescale. During the initial, slow motion stage, existing gravitational 
waves do not need to be resolved accurately in time, hence the timestep can be increased to match the hydrodynamic 
CFL condition. Monitoring of the fluid velocities, and adaptation of the time-step during the dynamic phase, will 
correctly capture the fast moving phase and the associated waveforms. 

As discussed in section W, the inner boundary for the hydrodynamical computation is taken to be further from 
the horizon than for the wave equation. In the coordinate system employed here (Boyer-Lindquist coordinates) the 
numerical integration of hydrodynamic equations becomes problematic close to the horizon, as the Lorentz factor 
diverges there. We checked that our results are not affected to a noticeable degree by the location of the boundary, 
but it is clearly desirable to have a uniform integration domain for the two systems, especially in that sensitive and 
interesting region. A satisfactory numerical extension of the coordinates inside the horizon, for the hydrodynamical 
equations, has been reported recently |4^-|4^. It appears that the computation of emission from more general flows 
around, in particular, rotating holes, will benefit from a complete transcription of the coupled problem onto regular 
coordinate patches (and non-singular tetrads). 

Radiative boundary conditions at infinity are delicate to impose on a finite grid, even for the linearized curvature 
perturbations that we are considering. Approximate schemes always achieve a partial transmission of wave energy 
across the boundary but this transmission can drop significantly for curved wave-fronts in a three-dimensional setting. 
Furthermore, imposing radiative conditions at a finite distance implies an effective truncation of the equation’s coeffi¬ 
cients. This interferes with physical features of the evolution which depend on terms beyond the principal part, such 
as quasi-normal ringing and tails. A developing programme for solving such problems is the “Cauchy-characteristic 
matching”. In this concept has been applied to equations arising in the context of black hole perturbation theory. 
For a review of this entire effort see . In the present work the problems with radiative conditions are exacerbated 
in an unexpected way: whereas for scalar wave equations incoming and outgoing radiation scales with the same power 
of r, for the spinorial equations considered here, the well known peeling properties of the Weyl tensor ensure that 
this is not so. Small amounts of reflection in outer boundaries, propagate back into the grid and get amplified as 
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they propagate inwards, entirely due to the radial scaling properties of incoming waves. In the absence of a seamless 
connection to infinity, an expensive, but safe, way to bypass reflection problems is to consider only the domain of 
dependence of the bounded initial hypersurface. 

The current implementation of the method allows the study of arbitrary axisymmetric flows onto a non-rotating 
hole. Using this setup we plan to extend the computations performed here to the case of matter motions with angular 
momentum and significant pressure support, as for example the dynamic collapse of the toroidal debris of a disrupted 
neutron star onto the black hole companion. Gravitational radiation emitted in this process would be the swan song 
of the once powerful binary emitter. With the parallel development of non-linear coupled codes (see e.g., p^ ), issues 
lying at the interface between perturbation theory and full numerical relativity become amenable to exploration, most 
interestingly the effects of self gravity. 
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FIG. 1. Calibration of the numerical error in computing the scattering and propagation of gravitational waves. Plotted is 
the percent relative error in the total emitted energy associated with a pure disturbance (i.e., no matter sources), after 400M 
of evolution. The insert illustrates the Cauchy converge to a reference grid of 800 x 20 points. The dotted line shows for 
comparison the theoretical second order convergence rate. The measured rate is about 2.2. The wave propagation converges 
relatively fast: fewer than 400 radial points within SOM from the black hole, are sufficient for a relative error of less than one 
percent. 
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FIG. 2. Calibration of the numerical error in computing wave generation by accreting matter. Plotted is the percent relative 
error in the total emitted energy associated with an infalling shell (see text for precise parameters) after 400M of evolution. 
The insert again focuses on the rate of Cauchy convergence of the total energy, computed with a reference grid of 2000 x 20 
points. The rate is measured to be about 1.9. The resolution requirements for achieving levels of error comparable to those 
of Fig. ^ are now considerable larger: 1200 zones inside 50M are needed to achieve an energy estimate accurate to about 2.5 
percent. 
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FIG. 3. Snapshots of density for a typical shell accretion computation. The symmetry axis 0 = 0 is horizontal, and the 
location of the black hole is depicted with the small black disk at the center. The initial, radially very compact profile of the 
density is shown at t = 0, followed by the gradual acceleration and infall shown at times t — lOOM and t — 200M. The 
final state is a relic (non-radiative) spherical pattern, corresponding to the low background density that “fills the vacuum” 
surrounding the shell. The t = 200M panel is a four-fold zoom into the interesting region. The lines are iso-density contours 
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FIG. 4. The gravitational wave time signal of a typical shell accretion event. The upper panel shows the quadrupole part 
of the wave signal normalized by the total mass of the shell, reaching a remote observer (r* = lOOM). The presence of two 
bursts is the first striking observation. The second burst is more than two times stronger, but the wavelength appears to be 
comparable. To bring out the full range of features, we show in the lower panel the logarithm of the same signal. The nature of 
the signal is now clearer: besides initial low frequency features ahead of each burst, both bursts are mostly in the quadrupole 
ringing frequency of the black hole (17M). Notice the transition region in between the two bursts; a non-oscillatory precursor 
wave is seen to rise steadily in amplitude, finally giving way to the strong oscillatory emission 
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FIG. 5. Spacetime diagram of wave and density evolutions. Light rays in the diagram are propagating at 45° angles. The 
domain depicted extends only to r* = 50M, in order to bring out the interesting inner region features. The field variables are 
sampled at the equator {9 = t:/2). Both panels show the contour lines of equal amplitude, for the wave signal (the logarithm of 
Z) on the left, and the logarithm of density on the right. In the right panel we see quite clearly the parabolic infall trajectory of 
the shell. The left panel shows the prodigious response of the black hole curvature to the presence of the shell. There are only 
three contour levels illustrated in the left panel, at logarithmic amplitudes of —6, 0, 2.5 respectively. The, lowest, —6 level is 
seen to form a light cone emanating from the initial location of the shell. The cone of influence reaches inwards to the horizon 
region, and excites ringing. The features of this ringing are traced summarily by the 0 contour level, e.g., the ringing is fully 
developed at around time t = TOM. As this initial burst subsides exponentially, the shape of the infalling shell is emerging, 
through the imprint it leaves on the curvature, as can be seen in the interval from 120M up until 210M. At t = 210M the 
shell accretes rapidly, and it quasi-static curvature pattern gives way to strong emission (traced by the 2.5 contour), which is 
seen to be here predominantly in QNM form. 
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FIG. 6. Spatial snapshots of wave profiles for a typical accretion event (part 1). The panels can be thought of as horizontal 
slices of the left column in Fig. |^, restoring one suppressed dimension at the corresponding time levels. The initial burst and 
associated ringing are shown at t = SOM and t = lOOM. At t = 150M those initial features have subsided sufficiently to reveal 
the quasi-static field surrounding the shell and being dragged with it towards the hole. At t = 200M, the density is near the 
peak of the potential (compare to Fig. ^), and the emission is reaching its climax. 
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FIG. 7. Spatial snapshots of wave profiles for a typical accretion event (part 2). At t = 250M the excitement is already 
gone: the remaining evolution through t = 400M only obtains the slow decay of the black hole quasi-normal mode ringing. 
The black hole returns to its quiescent state. 
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FIG. 8. Emitted energy as a function of the radial width of the infalling shell. More precisely, the energy dependence is 
shown with respect to the e-folding (proper) radial length for the shell density. All runs are for shells centered around r, = 35M, 
initially at rest, with the k parameter ranging from 0.1 to 100. The increased efficiency as a function of compactness is evident, 
and spans several orders of magnitude. In the extended limit the energy is shutting off, apparently as a power law of L* (see 
insert). The slope of the power law is htted well (dashed line) by a value of —2.4. In the thin shell limit, the energy asymptotes 
to a finite value which is about a third of the point particle limit. 
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FIG. 9. Emitted energy as a function of the initial velocity of the infalling shell. Shown is the maximum signal amplitude 
as a function of initial velocity, for velocities up to 0.1 times the speed of light. The infall velocity is seen here to have an 
important effect on the emission efficiency,; the amplitude for the largest initial velocity is more than twice than starting from 
rest. 
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FIG. 10. Excitations of the fundamental (1=2) black hole quasi-normal mode, as a function of the infalling shell width. 
The left column depicts the waveform for four select values of shell width (increasingly thinner shells from top to bottom). 
The large increase of per unit mass emission is evident. The wavelength of the emission is also increasingly shorter. The right 
column shows the logarithm of the signal, and brings out the qualitative changes emerging in the thin shell limit. More compact 
shells are seen to be first marginally (second panel), and then more clearly, exciting the QNM ringing. In the last panel the 
ringing is QNM frequency is dominant from t = 315M onwards, but it appears that the maximum amplitude is still emitted at 
a somewhat lower frequency. 
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FIG. 11. The density evolution for a shell collapse onto a moving black hole. The hole rushes up the symmetry axis (right 
side of the panels). By t = lOOM the initially spherical shell has been visibly distorted, as parts of it experience increased 
acceleration due to the proximity to the approaching hole. By t = 200M, the hole has gone right through the shell and we can 
see concentrated accretion to occur at the equatorial plane. At this time most of the gravitational wave emission has occurred. 
In the final snapshot (at t = 300M) we can see the slow accretion of the remaining shell material from the back side of the 
black hole. The lower panes progressively focus on the inner region. Despite the already astrophysically unlikely kick velocity 
of 0.1c, the energy released in this event (as I = 2 radiation) is measured to be only 2.9 x 10~^/M. The waveform is very 
similar to the bottom panel in Fig. 
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FIG. 12. The density evolution tor anisotropic, pancake , collapse. The shell collapses onto the black hole with the 
equatorial region falling in hrst, while the polar regions accelerate much later. By t = 150M the equatorial zones have accreted; 
further material keeps falling in from larger latitudes. The energy released here in the quadrupole mode is measured at 
2.4 X 10“®/r^/M, and comes predominantly from the equatorial accretion event. 










